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ABSTRACT 

We have studied the Bekenstein-Sandvik-Barrow-Magueijo (BSBM) model for the spatial and tem- 
poral variations of the fine structure constant, a, with the aid of full A-body simulations which 
explicitly and self-consistcntly solve for the scalar field driving the a-evolution . We focus on the 
scalar field (or equivalently a) inside the dark matter halos and find that the profile of the scalar field 
is essentially independent of the BSBM model parameter. This means that given the density profile of 
an isolated halo and the background value of the scalar field, we can accurately determine the scalar 
field perturbation in that halo. Wc also derive an analytic expression for the scalar-field perturbation 
using the Navarro-Frenk- White halo profile, and show that it agrees well with numerical results, at 
least for isolated halos; for non-isolated halos this prediction differs from numerical result by a (nearly) 
constant offset which depends on the environment of the halo. 



1. INTRODUCTION 

During the past ten years there has been continual in- 
terest in the possible time and space variation of fun- 
damental constants of Nature. This interest was stimu- 
lated by observations of quasar absorption spectra that 
are consistent with a slow increase of the fine structure 
"constant", a, over cosmological time scale (Webb et al. 
1999, 2001). Experimental and observational efforts to 
constrain the level of any possible time variation in fun- 
damental constants have a history that pre-dates mod- 
ern theories about how they might vary (for overviews 
see Barrow (2002, 2005), Uzan (2003) and Olive & Qian 
(2004)). Until quite recently all the observational stud- 
ies found no evidence for any variations. However, high- 
quality data from a number of astronomical observations 
have provided evidence that at least two of these con- 
stants, the fine structure constant: a = e 2 / he, and 
the electron-proton mass ratio: \i = m e /m p might have 
varied slightly over cosmological time. Using a data set 
of 128 KECK-HIRES quasar absorption systems at red- 
shifts 0.5 < z < 3, and a new many-multiplet (MM) 
analysis of the line separations between many pairs of 
atomic species possessing relativistic corrections to their 
fine structure, Webb et al. (1999, 2001) found the ob- 
served absorption spectra to be consistent with a shift in 
the value of fine structure constant, a, between those 
rcdshifts and the present day, with Aa/a = a(z) — 
a(0)/a(0) = -0.57 ±0.10 x 10~ 5 . A smaller study of 23 
VLT-UVES absorption systems between 0.4 < z < 2.3 by 
Chand et al. (2004) and Siranand et al. (2004) initially 
found Aa/a = — 0.6±0.6x 10 -6 by using an approximate 
version of the full MM technique. However, the reanaly- 
sis of the same data set by Murphy, Webb & Flambaum 
(2007, 2008) using the full and unbiased MM method 
increased the uncertainties and suggested a revised fig- 
ure of Aa/a = —0.64 ± 0.36 x 10 -5 for the same data. 
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These investigations relied on the statistical gain from 
large samples of quasar absorption spectra. Most re- 
cently, this observational programme has been extended 
to both hemispheres of the sky using KECK and VLT sam- 
ples of 153 absorption systems by Webb et al. (2010) 
and finds evidence consistent with an increase in a the 
northern sky but consistent with a slow decrease in a 
with time in the south. When combined these overlap- 
ping data sets are well fitted by a dipole with Aa /«o — 
(1.10±0.25) x 10 _6 r cos#, at measurement position r (rel- 
ative to Earth at r = 0) where 6 is the angle between the 
measurement and the axis of the dipole. These observa- 
tions suggest that we should develop an understanding 
of the spatial as well as the temporal consequences of 
varying constants. 

By contrast to these statistical searches for varying a, 
probes of the electron-proton mass ration can use sin- 
gle objects effectively. Reinhold et al. (2006) have found 
a 3.5<7 indication of a variation in the electron-proton 
mass ratio fi = m e /m pr over the last 12Gyrs: Afi/fi = 
(—24.4 ± 5.9) x 10~ 6 from H 2 absorption in a different 
object at z = 2.8. However, Murphy et al. (2008) have 
exploited the \i sensitivity of ammonia inversion transi- 
tions (Flambaum & Kozlov 2007) compared to rotational 
transitions of CO, HCN, and HCO + in the direction of 
the quasar B0218+357 at z = 0.68466 to yield a result 
that is consistent with no variation in fi when systematic 
errors are more fully accounted for: Afi/fi = (+0.74 ± 
0.47 sta t±0.76 syst cm) x 10~ 6 , corresponding to a time vari- 
ation of fi/ \i = (-1.2 ± 0.8 stat ± 1.2 By8tom ) x lO-^yr" 1 
in the best-fit ACDM cosmology. 

Any variation of a and fi today could also be con- 
strained by direct laboratory searches. These are per- 
formed by comparing clocks based on different atomic 
frequency standards over a period of months or years. 
Until very recently, the most stringent constraints on 
the temporal variation in a arose by combining mea- 
surements of the frequencies of Sr (Blatt et al. 2008), 
Hg+ (Fortier et al. 2006), Yb+ (Pcik et al. 2004), and 
H (Fischer et. al. 2004) relative to Caesium: a/ a — 
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(-3.3 ± 3.0) x lO-^yr" 1 . Cingoz et al. (2007) also re- 
cently reported a less stringent limit of a/a = — (2.7 ± 
2.6) x 10 -15 yr _1 ; however, if the systematics can be fully 
understood, an ultimate sensitivity of 10 _18 yr _1 is pos- 
sible with their method Nguyen et al. (2004). If a lin- 
ear variation in a is assumed with cosmic time then the 
Webb et al. (1999, 2001) quasar measurements equate to 
a/a = (6.4±1.4) x 10 _16 yr _1 . If the variation is due to a 
light scalar field described by a theory like that of Bckcn- 
stein, Sandvik, Barrow and Magucijo (BSBM from here on) 
(Bekenstein 1982; Sandvik, Barrow & Magueijo 2002), 
then the rate of change in the constants is exponen- 
tially damped during the recent dark-energy-dominatcd 
era of accelerated expansion, and one typically predicts 
a/a = 1.1 ± 0.3 x 10 _16 yr _1 from the Murphy et al. 
data, which is not ruled out by the atomic clock con- 
straints mentioned above. For comparison, the Oklo nat- 
ural reactor constraints, which reflect the need for the 
Sm 149 + n — > Sm 150 + 7 neutron capture resonance at 
97.3 meV to have been present 1.8 — 2Gyr (z = 0.15) 
ago, as first pointed out by Shlyakhter (1976), are cur- 
rently (Fujii et al. 2000) Aa/a = (-0.8 ± 1.0) x 10" 8 or 
(8.8±0.7) x 10 -8 (because of the double- valued character 
of the neutron capture cross-section with reactor temper- 
ature) and (Lamoureaux 2004) Aa/a > 4.5 x 10 -8 (6a) 
when the non-thermal neutron spectrum is taken into 
account. However, there remain significant environmen- 
tal uncertainties regarding the reactor's early history and 
the deductions of bounds on constants. The quoted Oklo 
constraints on a apply only when all other constants are 
held to be fixed. If the quark masses to vary relative to 
the QCD scale, the ability of Oklo to constrain variations 
in a is greatly reduced (Flambaum 2007). 

Recently, Rosenband et al Rosenband et al. (2008) 
measured the ratio of aluminium and mercury single- 
ion optical clock frequencies, /ai+/ fng+, over a period of 
about a year. From these measurements, the linear rate 
of change in this ratio was found to be (—5.3 ± 7.9) x 
10~ 17 yr -1 . These measurements provide the strongest 
limit yet on any temporal drift in the value of a: a/a = 
(—1.6 ± 2.3) x 10 _17 yr _1 . This limit is strong enough 
to conflict with theoretical explanations of the change in 
a reported by Webb et al. (1999, 2001) in terms of the 
slow variation of an effectively massless scalar field, even 
allowing for the damping by cosmological acceleration, 
unless there is a significant new physical effect that slows 
the locally observed effects of changing a on cosmological 
scales. Also, one expects inhomogeneous changes in a in 
scenarios where variations in a are induced by a heavy 
scalar field with a mass (m^). One would expect varia- 
tions in a on cosmological scales to differ from those on 
scales below the field's Jeans length, which is Oil/m^) 
(for a detailed analysis of global-local coupling of vari- 
ations in constants sec Clifton, Mota & Barrow (2005); 
Mota & Barrow (2004a,b); Shaw & Barrow (2006c,a,b); 
Olive & Pospelov (2008)). 

It has also been noticed that if 'constants' such as a 
or fi could vary, then in addition to a slow temporal 
drift one would also expect to see an annual modulation 
in their values. In many varying constant theories, the 
Sun perturbs the values of the constants by a factor 
roughly proportional to the Sun's Newtonian gravita- 
tional potential (Magueijo, Barrow & Sandvik 2002; 



Barrow, Sandvik & Magucijo 2002; Barrow & Mota 
2002) (the contribution from the Earth's gravitational 
potential is about 14 times smaller than that of the Sun's 
at the Earth's surface). Hence the 'constants' depend 
on the distance from the Sun. Since the Earth's orbit 
around the Sun has a small ellipticity, the distance, r, 
between the Earth and Sun fluctuates annually, reaching 
a maximum at aphelion around the beginning of July 
and a minimum at perihelion in early January. It was 
shown that in many varying constant models, the values 
of the constants measured here on Earth oscillate in a 
similar seasonal manner. Moreover, in many cases, this 
seasonal fluctuation is predicted to dominate over any 
linear temporal drift (Barrow & Shaw 2007). 

In this paper we will study the BSBM model for the 
spatial and temporal variations of the fine structure con- 
stant a, with the aid of full A-body simulations which 
explicitly and self-consistently solve for the scalar field 
driving a-evolution. We focus on the trend of the scalar 
field (or equivalently, a) inside the dark-matter halos, 
and find that the profile of the scalar-field fluctuation 
is essentially independent of the BSBM model parameter. 
This means that given the density profile of an isolated 
halo and the background value of the scalar field, we 
can accurately determine the scalar field perturbation in 
that halo. We also derive an analytic expression for the 
scalar-field perturbation using the Navarro-Frenk- White 
(NFW) halo profile, and show that it agrees well with 
numerical results, at least for isolated halos. For non- 
isolated halos this exact prediction differs from numeri- 
cal result by a (nearly) constant offset, depending on the 
environment of this halo. 

A brief outline of the remaining of this paper is as fol- 
lows: in § 2 we list the minimal set of necessary equations 
to understand the physics, and briefly describe our algo- 
rithm; § 3 displays the main numerical results and we 
then discuss and conclude in § 4. 

2. EQUATIONS AND ANALYSIS 

This section lists the equations which will be used in 
the iV-body simulations for the BSBM varying-a model 
Barrow & Mota (2003); Barrow, Magueijo & Sandvik 
(2002); Barrow, Sandvik & Magueijo (2002); 

Sandvik, Barrow & Magueijo (2002); Mota & Barrow 
(2004b,a). 

2.1. The Basic Equations 

The Lagrangian density for the BSBM model could be 
written as 



Cr. 



EM 



C r ,(1) 



where R is the Ricci scalar, k = 8ttG with G being the 
gravitational constant, ip is the scalar field; L m , £eMi C t 
represent respectively the Lagrangian densities for dust, 
electromagnetic field (including photons) and other radi- 
ation (such as neutrinos). The coupling function between 
the scalar field and the electromagnetic field in the BSBM 
model is e ^ 2 ^-v where ^/k is added so that y/~K(p = ip is 
dimensionlcss. In the simplest version of the model there 
is no potential for the scalar field. 
The dust Lagrangian for a point particle with mass mo 
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is 



ra 



£m(y) = ^^(y - x ) J gabioio 



(2) 



where y is the general coordinate and Xo is the coordinate 
of the centre of the particle. From this equation we derive 
the corresponding energy-momentum tensor: 



rpab 



m 



5(y-xo)±g±o- 



(3) 



Also, because ga^i^x^ = g a bU a u b — 1, in which u a is the 
four- velocity of the dark-matter particle centred at Xq, 
the Lagrangian can be rewritten as 



£ m (y) : 



mo 



<*(y - x ) 



(4) 



This result will be used below. 

Eq. (3) is just the energy-momentum tensor for a single 
matter particle. For a fluid consisting of many particles, 
the energy-momentum tensor will be 



rpab 




-Pcdmu u , 



(5) 

where V denotes a volume which is microscopically large 
but macroscopically small, and we have extended the 3- 
dimensional S function to a 4-dimensional one by adding 
a time component. Here, u a is the averaged four-velocity 
of the matter fluid. 
Using 



rpab 



5g a b 



(6) 



it is straightforward to show that the energy-momentum 
tensor for the scalar field is 



rjVf> 



ab 



(7) 



Therefore the total energy-momentum tensor is 
Tab = V a (pV b (y5 - ig a &V c </?V V 



t: 



EM 



(8) 



where T™£ = p m u a Ub, T* b is the energy-momentum ten- 
sor for radiation fields except photons, and T® b M for pho- 
tons. The Einstein equations are 



G a b — K T a 



(9) 



in which G a b = R a b — \gabR is the Einstein tensor. Note 
that due to the extra coupling between the scalar field, 
tp, and the electromagnetic field, the energy-momentum 
tensors for either will no longer be separately conserved, 
but instead we have 

V&T^ = (g ab C E M + ?em) V b <p. (10) 

However, the total energy-momentum tensor is certainly 
conserved. 

Meanwhile, the scalar field equation of motion is 



EM. 



(11) 



where □ = V a V a . This equation governs the time evo- 
lution and spatial configuration of the scalar field. 

Eqs. (8, 9, 10, 11) summarize all the physics needed 
for the following analysis. However, when making use of 
them we should also specify the form of the electromag- 
netic matter. For example, if it is a pure electromagnetic 
field (photons), then we have £em = \ (E 2 - B 2 ) = 
in which E, B stand for the electric and magnetic fields. 
Thus from the time component of Eq. (10) we obtain the 
(background) evolution equation for photon density as 



AHp r = 2ipp r 



(12) 



where remember that tp = y/nip. 

It then seems that, by the same reason, the right 
hand side of Eq. (11) also vanishes, leaving the scalar 
field unsourced. This may not be true however, as non- 
relativistic matter could also contribute to Cem and thus 
T^ b M . For example, in baryonic matter £em ~ E 2 /2, and 
for neutrons and protons this electromagnetic contribu- 
tion to the total mass can be of order 10 -4 ; in supercon- 
ducting cosmic strings Cem ~ —B 2 /2 where pem ~ B 2 /2 
so that |£em/pem| ~ 1- In the BSBM model, in order to 
simplify the situation, it is assumed that Cem/ Pm = C 
where £ is a constant with a modulus between and 
~ 1, either positive or negative, and p m is the density for 
non-relativistic matter. 

Thus the scalar field equation gets sourced by a term 
proportional to 

□<^=-2V7<e- 2 ^Vm- (13) 

Since the part of £em which affects the scalar field is a 
constituent of the non-relativistic matter and is presum- 
ably moving with the matter particles, we could combine 
Eq. (10) and the conservation equation for the dust mat- 
ter (no including electromagnetic contribution) to write 
a new conservation equation for the particle: 

V b T« fe +EM = 2^ {g ab C EM + Til) V 6 V- (14) 

Although we have assumed above that £em = (p m , we 
still lack a knowledge about Tg^ , whose relation to Cem 
could be complicated. Here for simplicity we assume that 
Tem = — (PmU a u b . Then it is easy to find that the time 
component of this equation reads 

Pm+EM + 3Hp m+ EM = (15) 

while the i-th spatial component of it gives the following 
(modified) geodesic equation 



ab^O 



b a = 2Cv^ (g lb - u l u b ) 



(16) 



where xq is the coordinate of the centre of a parti- 
cle, and the right hand side represents a fifth force on 
the particle Li & Zhao (2009, 2010). The assumption 
Tem ~ ~(,PmU a u b might seem unappealing, but as we 
shall see below, because |£| <?C 1, the fifth force is much 
weaker than gravity and thus has negligible effects in the 
clustering of matter in any case; ultimately it is only the 
BSBM assumption £em = CPm that is important in theo- 
retical predictions of the spatial and temporal variations 

of a, given that a — e 



2.2. Analytical Approximation 
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The scalar field equation of motion, which controls the 
dynamics of the scalar field ip, is generally complicated, 
because it depends nonlincarly on p, which both evolves 
in time and fluctuates in space. Fortunately, for the ma- 
jority of applications the scalar field potential and cou- 
pling function are not nonlinear enough to give the scalar 
field a very heavy mass so as to make it fluctuate strongly. 
The nice thing about this is that in certain places of 
the equations one may then forget the scalar field per- 
turbation and simplify these equations accordingly. In 
Li & Barrow (2010), for example, it is shown that such 
simplification is very good an approximation (see how- 
ever Li & Zhao (2009, 2010) for an opposite extreme for 
which the scalar field potential is very nonlinear so that 
such simplification does not work). 

In the BSBM model, there is no scalar field potential 
and the coupling function is close to linear if V^M ^ 1 
(which is the case for our interested parameter space 
(Barrow, Magueijo & Sandvik 2002)). We therefore ex- 
pect the fluctuation of ip to be very weak and assume 
that y^I'VI *C v^M <^ 1 (which we shall confirm below 
using numerical simulations). Under such an assumption 
the Poisson equation could be written as 



V 2 $ = 47rGa 3 ^ 



-4:irGa 3 p ri 



1 + Ce 
14- (e 



•AirGa 3 (pm - p m ) 



(17) 



where p the background value of p and Sp its perturba- 
tion; PmiPm are respectively the local and background 
matter density; $ is the gravitational potential and a 
the cosmic scale factor; V x is the derivative with respect 
to the comoving coordinate x. To obtain Eq. (17) we have 
used the fact that in the BSBM model \C,e~ 2 ^\ < 1. Note 
that Eq. (17) clearly indicates that the gravitational po- 
tential is essentially not influenced by the scalar field p. 

For the scalar field equation of motion Eq. (13), be- 
cause the background part (which has no spatial depen- 
dence) can be solved easily, so we subtract that from the 
full equation to obtain an equation of motion for 5p only 
(remember that p = p + Sep). Furthermore, we drop the 
time derivative terms of Sp> as they are small compared 
with the spatial gradients (i.e., work in the quasi-static 
limit). The final equation for Stp then becomes 



V 2 (ay/ndp) = 2(k 



-2-\/k<p 



^16irGa 3 (( Pm - p m ) (18) 

where we have used k = 8irG and ^/k|<5<^| <C \/^M -C 1. 

Comparing Eqs. (17, 18), it is evident that the source 
terms are the same up to a constant coefficient 4£. Con- 
sequently, we shall have 

av ^^(x)«4C$(x). (19) 

Note that this equation, together with the geodesic equa- 
tion Eq. (16), implies that the magnitude of the fifth force 
(force due to exchange of scalar field quanta between par- 
ticles) |f| satisfies 



C|V(av^MI ~ 4C 2 |V$| 



(20) 



Therefore the ratio between the magnitudes of the fifth 
force and gravity is of order ( 2 < 10~ 12 — 10~ 8 <C 1. This 



Fig. 1. — The time evolution of ip for the four models considered 
in this work, with ( = -2x (solid curve), 2 X 1CT 6 (dotted 

curve), — 5x 10 — 6 (dashed curve) and 5x 10 — 6 (dash-dotted curve) 
respectively. The horizontal axis is the cosmic scale factor a(t) and 
the vertical axis plots y/n(p in unit of 10~ 4 . 

implies that the fifth force cannot influence the structure 
formation significantly 

3. SIMULATIONS AND RESULTS 
3.1. p Perturbation vs. Gravitational Potential 

To study in details the behaviour of the scalar field and 
hence the fine structure constant a in the BSBM model, 
we have performed A-body simulations for four different 
models, with ( = ±2 x 10 -6 and ±5 x 10 -6 respectively. 
The physical parameters we adopt in all simulations are 
as follows: the present-day dark-energy fractional energy 
density 17de = 0.743 and Q m = Ocdm 4- = 0.257, 
H = 71.9 km/s/Mpc, n s = 0.963, er 8 = 0.761. These are 
in accordance with the concordance cosmological model 
preferred by current data sets. Our simulation box has a 
size of 64ft," 1 Mpc, in which h = ff /(100 km/s/Mpc). 
In all those simulations, the mass resolution is 1.114 x 
10 9 /i _1 Mq; the particle number is 256 3 ; the domain 
grid {i.e., the coarsest grid which covers the whole simu- 
lation box) has 128 3 equal-sized cubic cells and the finest 
refined grids have 16384 cells on each side, corresponding 
to a force resolution of ~ 12/i _1 kpc. Detailed descrip- 
tion about the A-body simulation technique and code 
for the (coupled) scalar field models could be found in 
Li & Zhao (2009, 2010) and will not be presented here. 

Because the BSBM model (with the parameter £ con- 
strained by data) involves a very weak coupling between 
matter and the scalar field p, the presence of the latter 
and its coupling have negligible influences on the back- 
ground (ACDM) cosmology, although the opposite is not 
true since there is no scalar field potential and thus the 
dynamics of ip is controlled entirely by the coupling. In 
our simulations, we compute the full background cos- 
mology and evolution of ip on a predefined time grid us- 
ing MAPLE, and interpolate to obtain the corresponding 
quantities which are needed in A-body simulations. De- 
tails of this procedure can be found in the Appendix C 
of Li & Barrow (2010), and Fig. 1 shows the background 
evolution of y/lvp for the 4 models considered here. Ob- 
viously the condition ^JUp <C 1 is satisfied, justifying the 
approximation we used above to derive Eq. (19). 
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Fig. 2. — (Color Online) Scatter plot of the scalar field perturbation, ay/nSip, (vertical axis, in unit of 10 — 11 ) vesus gravitational potential 
<t> (horizontal axis) in code units. The black solid line is the analytical approximation adip oc <I>; each green dot represents the corresponding 
result measured in a cell from a particular slice that is randomly selected from the simulation box. The columns are for four models with 
£ = ±2, ±5 X 10 — B as shown above each panel, and the rows are for three different output times a = 0.3,0.5, 1.0, also shown above each 
panel. Note that the slopes for the solid lines differ because of the different values of £ . 



One nice thing about the BSBM model is its sim- 
plicity, and it turns out that the background evolu- 
tion of ip (and thus a) in different cosmic epochs 
can be well described by some analytical formulae 
Barrow, Magueijo & Sandvik (2002). Therefore in the 
present study we shall mainly focus on the spatial varia- 
tion of ip and a (especially in virialized halos) . 

As mentioned above, because there is no potential for 
the scalar field ip and because \fnip <C 1 for our choices 
of £, so the scalar field equation of motion (in the quasi- 
static limit) for Stp and the Poisson equation share the 
same source up to a constant coefficient, and therefore 
we expect a8<p ex $ across the whole space. This is what 
we have found in Li & Barrow (2010) for a different cou- 
pled scalar field model where the scalar field potential is 
negligible. Indeed, this could serve as a test of the scalar 
field solver in our TV-body simulation code. 

To check that our code does recover the analytical ap- 
proximation, we have plotted in Fig. 2 the comparison of 
the scalar field perturbation ay/RSip and gravitational po- 
tential $ from a slice of the simulation box. As indicated 
by this figure, the agreement between the numerical re- 
sults (green dots) and analytical approximation (black 
solid line) is remarkably good, implying that the numeri- 



cal code works well. Therefore to a high precision we can 
assume that a^/nSip oc $ everywhere, a fact which shall 
be used below to obtain an analytical expression of Sip 
in dark matter halos. 

3.2. Spatial variation of ip in Halos 

In the standard picture, the galaxies where observers 
live generally locate inside the dark matter halos, which 
to the simplest approximation are just spherical clusters 
of matter with a universal NFW (Navarro et al. 1996) 
density profile. 

We are certainly interested in the (possible) variation 
of a inside the halo we reside in. For example, there 
has been a great deal of analytical work on how signifi- 
cantly the local value of a could deviate from its cosmo- 
logical counterpart (Mota & Barrow 2004b, a; Jacobson 
1999; Wetterich 2003; Shaw & Barrow 2006c,a,b). Also, 
if the spatial variation of a is strong enough, then it 
might have impact on our observation of the spectra for 
the stars from our Galaxy and other galaxies. 

From our simulation output, it is easy to identify the 
dark matter halos (Knebe & Gibson 2004; Li & Barrow 
2010). Now what we want to do here is to measure the 
quantity yfkSip as a function of distance R to the halo 
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Fig. 3. — The profile of a^/Ti&ip inside the halos for the four models with C, = ±2, ±5 X 10 as indicated at the top of each panel. 
We have considered the 80 most massive halos from the simulation box, and divided them into six bins with M/ (l0 14 M Q ) 6 [2.0, oo), 
[1.0, 2.0], [0.6, 1.0], [0.3, 0.6], [0.2, 0.3] and [0.1, 0.2] respectively. Each curve represents the averaged profile of a^fKhip for one bin (the more 
massive bins always correspond to larger, either positive or negative - depending on the sign of f — deviations from zero). All results are for 
output time a = 1.0 (today). The horizontal axis is the distance R to the halo centre, in units of h~ kpc and the vertical axis is a^/nSip 
in units of 10~ 10 . 



centres, assuming that the halos are exactly spherical. 
For this we have recorded the value for ^JH8tp at the po- 
sition of each particle, and then presumably we could 
divide each halo into a number of spherical shells, deter- 
mine the radius of each shell and compute the average 
value of yfnSip in the shells. However, there is some sub- 
tlety in the computation of the averaged y/nSip. 

The problem is that, since we have only recorded the 
information for ^/nSip at the positions of the particles, we 
do not have a fair sampling points. Because y/HSip tends 
to be different in regions of different particle number den- 
sities, the high density regions will be over-sampled and 
low density regions under-sampled, resulting in a bias as 
we are trying to determine the spatially averaged, rather 
than the particle averaged value, of i/nSip. To test how 
big the bias could be, we use an approximation as follows: 
firstly we divide each spherical shell into N equal-sized 
volumes which are small enough so that the particle num- 



ber density do not change much inside each of them, then 
we compute the particle average of t/k5(p in each of these 
volumes, call it (^/KS(p)i for i = 1 ■ ■ • N, then the spatial 
(volume) average of y/nSip is given by 

1 N 

(^)vol« ^2(y/K,Sip)i. (21) 

i=l 

More precise treatments of the volume average could be 
obtained by using space tessellations, such as Delaunay 
triangulation, but this is too technical and thus beyond 
the scope of this work. Anyway, using our approximation 
Eq. (21), we find that the bias caused by using particle- 
number average is at most 1 ~ 2%, which is not unac- 
ceptable considering that the sphericity of halos is al- 
ready an approximation. 

Fig. 3 shows the profile of ^/kSlp inside the dark matter 
halos. Instead of plotting this halo by halo, we have se- 
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(=-2e-6, a=1.0 (=2e-6, o=1.0 
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Fig. 4. — As Fig. 3, except that this figure shows the profiles of Sip/<p instead of a^fnhip. 



lected the 80 most massive halos from our simulation box, 
divide them into six bins with Mj (10 14 Mq) G [2.0, oo), 
[1.0,2.0], [0.6,1.0], [0.3,0.6], [0.2,0.3] and [0.1, 0.2]. Then 
we compute the averaged profile of yfJiSip in each of the 
six bins (the six curves in Fig. 3). As expected, the larger 
the halo is, the deeper the gravitational potential is and, 
because a^JHSpi « <f>, the larger Iv^^l is. Meanwhile, go- 
ing from the halo centre outwards one finds that 
gradually decreases (towards 0), which makes sense be- 
cause yfnSip — >• as R — s- oo. Finally, we also notice that 
the value of ^JUSip depends sensitively on Q and |-\/K<5y>| 
increases with |£|. 

Fig. 3 also shows that for our choices of £ the value of 
\^/k6<p\ is typically of order 10 -10 ~ 10~ 9 , which is quite 
small but gives us no idea how large the fluctuation in 
ip is. For the latter we have instead plotted the profile 
of the quantity Sip/ if inside the 80 halos distributed in 
six bins. Interestingly, unlike Sip, the quantity Sip/ip does 
not depend on the sign of (, and indeed it almost does 
not depend on the magnitude of £ cither! This, together 
with the facts that (a) ay/KSip(R,a) cx $(i?, a) and (b) 
the presence of the scalar field and its coupling to matter 



have negligible effect in the structure formation (so that 
the halo density profile remains NFW) , indicates that the 
fluctuation of ip in the BSBM model only depends on the 
non-BSBM physical parameters, and once we have solved 
the background value of tp we might get some pretty good 
idea about the ^fnSip profile without solving it explicitly 
(see below for further details). Back to the size of Sip/ip, 
from Fig. 4 we see clearly that it is of order 10 ~ 6 ~ 10 -5 , 
which is too tiny to produce any observable effects in the 
spatial variation of a. 

Now, given that the NFW profile for dark matter ha- 
los is (expected to be) preserved in the varying-a simu- 
lations, we wonder whether it is possible to derive some 
analytical (approximate) formula for the profile of ^fnip. 
For this let's recall that the NFW profile is expressed as 



where p c is the critical density for matter, /? is a di- 
mcnsionlcss fitting parameter and R s a second fitting 
parameter with length dimension. f3 and R s are gener- 
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Fig. 5. — The fitting curves for the circular velocity, V c , in the halo using the parameterization Eq. (23). We show the results for nine 
halos selected from the 80 most massive ones in the simulation box, and their masses are given near the bottom of the each panel (the 
mass is decreasing from the upper-left corner to the lower-right corner). Solid curves are the fitting formulae while crosses are the JV-body 
simulation results. All results are at the output time a = 1.0; the horizontal axis is the distance from the halo centre in unit of h~ 1 kpc 
and vertical axis denotes V c , in unit of km/s. Note that only the model with f = —2 X 10~ 6 is displayed for clarity but other models give 
similar results. 



ally different for different halos and should be fitted for 
individual halos. 

We have checked the halos in our analysis and found 
that the majority of them are indeed very well described 
by Eq. (22), confirming our earlier argument that the 
coupled scalar field effect is too tiny to change the struc- 
ture formation 1 . However, we shall not use the fitting to 
Eq. (22) in this work, mainly for two reasons: first, the 
dark matter density profile is in general difficult to mea- 
sure directly or precisely, while in contrast the circular 
velocities V c of the stars rotating about the halo centre 
are easier to measure; second, V c as a function of radius 
R is more closely related to t/k8ip(R), which will become 
clear later. As a result, we shall use fittings to V c from 
here on. 

Assuming Eq. (22) as the density profile and sphericity 
of halos, we could derive V c easily as 



V C 2 (R) = 



GM(R) 



R 

= A-kGPp c R? s 



R 



In 1 



R 

Ri, 



1 



R s + R 



(23) 



1 Indeed the NFW profile is quite robust, and even the scalar field 
does influence the structure formation significantly it is often still 
preserved. See examples for the coupled quintessence (Baldi et al. 
2010), RcBEL (Keselman, Nusser & Peebles 2010) and extended 
quintessence (Li, Mota & Barrow 2010) models. 



where M(R) is the mass enclosed in radius i?, and again 
this is parameterized by j3 and R s . From a simulation 
point of view, it is straightforward to measure M(R) and 
then fit (3 and R s ; from an observational viewpoint, it is 
easy to measure V C (R), which could again be used to fit 
(3 and R s . 

To show how good the fittings could be, we pick out 9 
halos with different masses and sizes from our simulation 
box, fit the corresponding and R s using the measured 
M(R), and plot these in Fig. 5. As can be seen there, 
the fitting results (solid curves) agree with the simulation 
results (crosses) quite well (in particular for the halo with 
M = 3.895 x 10 13 il/ o ). 

To see how this could be related to y/~Ktp, we remember 
that in the above we have shown that a^/Tup = 4£$, and 
so all we need to do is to find an expression for $(-R). For 
this, we use the fact that the potential inside a spherical 
halo is given as 



GM(r) 



dr + C 



(24) 



in which GM(r)/r 2 is the gravitational force and C is a 
constant to be fixed using the fact that <f>(i? = oo) = $00 
where i>oo is the value of the potential far from the halo. 
Using the formula for GM(r)/r 2 given in Eq. (23) it is 









Fig. 6. — (Color Online) The nalytic approximation compared with the numerical simulation results for the profile of a^JUbip in the same 
nine halos as in Fig. 5. Purple crosses are the numerical results, the solid curve represents the analytical approximation Eq. (29) with 
3?* = 'I'oo = 0, and the dashed curve denotes Eq. (29) with some tuned value of <E>„. The parameters f} and R s are the best-fit values in 
Fig. 5. All results are at the output time when a = 1.0; the horizontal axis is the distance from the halo centre, in units of h~ 1 kpc and 



the vertical axis denotes a^/nStp, in units of 10 
give similar results. 



Note that only the model with f = —2 X 10 is displayed for clarity but other models 



not difficult to find that 
nfl GM{r) 



-dr = ATiGfip c ir s 



lo i 
and so 

C = $ c 
Then it follows that 



f 

R~s 



4irGf3p c R 2 s 



$(i?) = $ c 



4^G/3p c ^ln| 



R 

R, 



If the halo is isolated, then <&oo = and we get 
$(i?H-47rG/3 Pc ^ln (l + f- 



(25) 
(26) 

(27) 



However, in iV-body simulations, we have a large number 
of dark matter halos and no halo is ideally isolated from 
the others. In such situations, $oo in Eq. (26) should be 
replaced by 



$(i? = .R* >i? vir )^0 



(28) 



where i?„ is some radius large compared with i? v i r (the 
virialized halo radius) but small compared with inter- 



halo distances. Then we get 



ay/K6(p(R)=4£ 



-47rG/?p c ^ln(\ + J- 



(29) 



As an example to show how well Eq. (2!)) works, we 
show in Fig. 6 the results of a\JTi8tp(R) for the same halos 
used to fit /3 and R s in Fig. 5. Here the crosses represent 
the values of a^JTiSip measured from the iV-body simu- 
lations and the curves our analytical approximations, in 
which the solid curve is obtained by setting $* = $oo = 
while the dashed curve is from tuning appropriately 
Obviously in most cases there is a (nearly) constant shift 
of the approximation with respect to the numerical re- 
sults, which accounts for the nonzero- ness of $*. 

Note that Eq. (29) captures the shapes for ay/ESip in 
various halos, but there is still one free parameter to 
be tuned to match the numerical results. This parameter 
summarizes our lack of knowledge about the environment 
which the considered halo resides in. As a result, the 
formula Eq. (29) is most suitable to apply in isolated 
halos, while for residential halos some extra work remains 
to be done to make it accurate. 

Alternatively, one could also consider Eq. (29) as a 
3-paramctcr parameterization of ay^Sip, for which the 
three parameters /3, R s and <&* could be fitted using the 
results from A^-body simulations. Given all the above re- 
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suits, we expect that this could produce some nice fitting 
curves too, but we shall not expand on this point here. 



4. SUMMARY AND CONCLUSION 

To summarise, in this paper we have studied the be- 
haviour of the BSBM varying-a model in the highly non- 
linear regime of large scale structure formation, with the 
aid of full A-body simulations that explicitly solves the 
scalar field which controls the temporal and spatial vari- 
ations of a. 

We have checked that, because of the weak coupling to 
matter and the lack of (nonlinear) potential, the scalar 
field is indeed very light everywhere and thus does not 
cluster significantly, i.e., the spatial fluctuation of ip is 
tiny. Because of this property, we have been able to sim- 
plify the field equations, which in turn suggest that the 
scalar field perturbation Sip is proportional to the gravi- 
tational potential $ [cf. Eq. (19)]. The numerical simula- 
tions then conform that such a simplification is justified 
to high precision. 

We then concentrate on the profiles of scalar field in- 
side virialized halos, which galaxies (and observers) are 
supposed to reside in. Figs. 3 and 4 display the averaged 
profiles of a^JHSip in the most massive halos from our 
simulation box, and they show that a^/nSip decreases as 
one goes from the halo centre outwards. In addition, the 
heavier the halo is, the deeper the gravitational potential 
will be and, as a result of Eq. (19), the larger ayfn\Sip\ 
is. Interestingly, although Sip does depends on the value 
of the BSBM parameter (, the quantity Sip/ip is essentially 
independent of it (i.e., the same for BSBM models with 
different £ from our interested parameter space). 

Thanks to the smallness of the scalar field coupling, the 
background expansion rate and the source to the Pois- 
son equation are essentially unaffected by the scalar field, 
while at the same time the fifth force is much weaker than 
(~ C 2 times) gravity so that its effect is also negligible. 
As a result, the structure formation itself is indistinguish- 
able from that in pure ACDM model. In particular, the 
halo density profiles are very well described by the NFW 
fitting formula. Furthermore, the circular velocities V c 
inside halos from the simulations are also well fitted us- 
ing the analytical formula for V c derived assuming NFW 
profiles [cf. Fig. 5]. 

As another check of the fact that the scalar field pertur- 
bation Sip is proportional to the gravitational potential 
$ [cf. Eq. (19)], we have derived an analytical expression 



for a^/nSp again by assuming NFW halo density profiles 
[cf. Eq. (29)]. We adopt the NFW parameters fitted using 
the circular velocities measured from simulation outputs 
in Eq. (29) and make predictions for a^fnSip in different 
halos. As shown in Fig. 6, the predictions agree very well 
with the numerical results for ay/HSip, if we take into ac- 
count the fact that halos are generally not isolated but 
living in potential wells produced by other halos. 

Our results suggest that for simple coupled scalar field 
models such as BSBM (and the one studied in Li & Barrow 
(2010)), the properties of the scalar field perturbation 
could be studied without solving the scalar field equa- 
tion of motion explicitly (which is time-consuming). We 
could cither extract them from say the halo density pro- 
files [using our Eq. (29)] of ACDM A-body simulations, 
or from the data on the galaxy rotation curves [using 
Eqs. (23, 29)] from observations. In the latter case, it is 
interesting that two seemingly uncorrelated things could 
be studied together and only once. 

Let us stress that the above nice properties are only 
present because of the smallness of the fluctuation in the 
scalar field, which is in turn due to the lack of significant 
nonlinearity in its equation of motion. If the scalar field is 
a chameleon, then its spatial fluctuation could be strong 
(Li & Zhao 2009, 2010; Zhao et al. 2010) and it then be- 
comes impossible to obtain a simple analytical formula 
for ay/tiSip, such as Eq. (29). A-body simulations will be 
the only tool to study such models, which we hope to 
investigate in the future. 
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